\chapter{Numerical Modeling of the All-Optical Switch}
\label{cha:nummodel}

In Chapter~\ref{cha:patterns}, I presented a theoretical model that describes pattern formation in counterpropagating beam systems. Chapters~\ref{cha:patterns_in_Rb} and \ref{cha:switch} describe my experimental observations of pattern formation and all-optical switching, respectively, in such a system where warm rubidium vapor serves as the nonlinear material. In parallel to my experimental work, I seek to determine the minimum ingredients required for a theoretical model to describe the switching behavior I observe experimentally. Previous attempts to simulate numerically all-optical pattern formation have been successful. Based on the model of Firth and Par\'e, numerical simulations performed by Chang \etal describe hexagonal pattern formation in a counterpropagating beam system \cite{Chang_1992aa}. This model, and the results of Chang \etal, serve as a starting point for my own numerical simulations.

The purpose of this Chapter is to describe my extension of \cite{Chang_1992aa} to simulate all-optical switching with transverse patterns. Specifically, I simulate the effect of a weak switch beam on the orientation of the hexagonal pattern generated by gaussian pump beams that counterpropagate through a Kerr-type nonlinear medium. Simulations of the time response of this system show behavior that is qualitatively similar to my experimental observations. In particular, the response time increases as the switch-beam power decreases.

In regards to the content of this Chapter, I owe extra thanks to Dr. Alex Gaeta of Cornell University for providing me with working Fortran code that implemented the split-step beam-propagation method in one transverse dimension (see Sec.~\ref{sec:split_step_beam_propagation}). Although none of the original code remains in my simulations, without getting started by standing on his shoulders, I would not have been as readily successful.

\section{3-D nonlinear model} % (fold)
\label{sec:3_d_kerr_model}
In order to simulate a full three-dimensional system, the model of Firth and Par\'e, presented in Sec.~\ref{sec:transverse_patterns}, is generalized here to include both transverse dimensions. This model is scalar, \emph{i.e.}\ it does not account for the vector nature of the fields, and hence cannot describe polarization instabilities, and it also does not include absorption effects. Nonetheless, it is sufficient to describe pattern formation in counterpropagating-beam nonlinear optical systems. As presented in Chapter~\ref{cha:patterns}, the following two equations describe the forward and backward fields counterpropagating through a Kerr medium, (also see Eq.~(\ref{eqn:field_amps})),
\begin{align}
		\left(\frac{\partial}{\partial z}+\frac{n_0}{c}\frac{\partial}{\partial t}\right)F&=\frac{i}{2k}\frac{\partial^2}{\partial x^2}F+i(|F|^2+2|B|^2)F\label{eqn:kerr_cprop_F}\\
		\left(-\frac{\partial}{\partial z}+\frac{n_0}{c}\frac{\partial}{\partial t}\right)B&=\frac{i}{2k}\frac{\partial^2}{\partial x^2}B+i(|B|^2+2|F|^2)B,\label{eqn:kerr_cprop_B}
\end{align}
where $n_0$ is the linear refractive index in the medium, $c$ is the speed of light in vacuum, $k$ is the wavevector within the medium, and $F$ ($B$) is the forward (backward) field amplitude. These equations further assume a nonlinear medium with positive $n_2$ as appropriate for my experimental system, described in Chapter~\ref{cha:switch}.

To normalize the variables in this model, I choose to scale the longitudinal units by the cell length $L$, the transverse units by the beam waist $w_0$, and hence the field strength carries additional length units as illustrated in Table.~\ref{table:normalized}, where $t_r=n_0 L/c$ is the transit time through the medium and the Fresnel number of the system ($\mathscr{F}=w_0^2 /\lambda L$) relates the transverse and longitudinal length scales. The Fresnel number can also be considered as a measure of the number of transverse modes supported by the geometry of the system. For large $\mathscr{F}$, many transverse modes are allowed and rich spatial structure is typically observed. The numerical results presented in this chapter are for $\mathscr{F}=5.3$, corresponding to the experimental results presented in Chapter~\ref{cha:switch}.

\begin{table}[htbp]
  \begin{center}
	\begin{tabular}{c|c}
	Normalized&Original\\
	\hline
	$\tilde{z}$&$z/L$\\
	$\tilde{x}$&$x/L$\\
	$\tilde{t}$&$t/t_r$\\
	$\tilde{F}$&$F\sqrt{L}$\\
	$\tilde{B}$&$B\sqrt{L}$\\
	$\kappa$&$Kw_0$\\
	\end{tabular}
  \end{center}
	\caption[Normalized units for Kerr model]{Equivalence between normalized units and the original physical units.}
	\label{table:normalized}
\end{table}

Rewriting Eqs.~(\ref{eqn:kerr_cprop_F}) and (\ref{eqn:kerr_cprop_B}) using these variables yields the following pair of dimensionless equations:
\begin{align}
		\left(\frac{\partial}{\partial \tilde{z}}+\frac{\partial}{\partial \tilde{t}}\right)\tilde{F}&=\frac{i}{4\pi\mathscr{F}}\frac{\partial^2}{\partial \tilde{x}^2}\tilde{F}+i(|\tilde{F}|^2+2|\tilde{B}|^2)\tilde{F}\\
		\left(-\frac{\partial}{\partial \tilde{z}}+\frac{\partial}{\partial \tilde{t}}\right)\tilde{B}&=\frac{i}{4\pi\mathscr{F}}\frac{\partial^2}{\partial \tilde{x}^2}\tilde{B}+i(|\tilde{B}|^2+2|\tilde{F}|^2)\tilde{B}.
\end{align}
Generalizing this model to two transverse dimensions is straightforward. The transverse derivative in $x$ is simply replaced by the transverse Laplacian operator 
\begin{equation}
  \nabla_\perp^2=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}.
\end{equation}
For simplicity, I will drop the tilde notation, and for the remainder of this Chapter, the variables are assumed to be in their normalized units, unless otherwise noted. The amplitude equations, including two transverse dimensions, are thus given by
\begin{align}
		\left(\frac{\partial}{\partial z}+\frac{\partial}{\partial t}\right)F&=\frac{i}{4\pi\mathscr{F}}\nabla_\perp^2F+i(|F|^2+2|B|^2)F\label{eqn:model_eqns_F}\\
		\left(-\frac{\partial}{\partial z}+\frac{\partial}{\partial t}\right)B&=\frac{i}{4\pi\mathscr{F}}\nabla_\perp^2B+i(|B|^2+2|F|^2)B.\label{eqn:model_eqns_B}
\end{align}
These amplitude equations can now be integrated by numerical methods. The following section describes a technique for reducing this set of PDEs to a pair of ODEs, and the following section describes a particularly efficient numerical integration scheme known as the split-step beam-propagation method (BPM) \cite{Fleck_1976aa, Moloney_1982aa, Burzler_1996aa}.

% section 3_d_kerr_model (end)

\section{Method of Characteristics} % (fold)
For a first-order hyperbolic partial differential equation (PDE), such as Eqs.~(\ref{eqn:model_eqns_F},\ref{eqn:model_eqns_B}), the method of characteristics can be used to find curves, known as characteristics, along which the behavior of the PDE is simplified and can be described by an ordinary differential equation (ODE). To solve for the characteristics of Eqs.~(\ref{eqn:model_eqns_F},\ref{eqn:model_eqns_B}), and for notational simplification that will become clear in later sections, it is useful to re-write the amplitude equations in terms of two operators $\mathscr{L}$ and $\mathscr{N}$
\begin{equation}\label{eqn:char_wave}
		\left(\frac{\partial}{\partial z}+\frac{\partial}{\partial t}\right)F=(\mathscr{L}+\mathscr{N})F,
\end{equation}
where 
\begin{align}
	\mathscr{L}(F)&=\frac{i}{4\pi\mathscr{F}}\nabla_\perp^2F,\label{eqn:linear_op}\\
	\mathscr{N}(F)&=i(|F|^2+2|B|^2)F,\label{eqn:nonlinear_op}
\end{align}
are operators describing linear diffraction and nonlinear wave mixing respectively, with similar equations for $B$. To apply the method of characteristics to Eq.~(\ref{eqn:model_eqns_F}), I start by parameterizing $z$ and $t$ as $z=z(s), t=t(s)$ for the parameter $s$ and assume an initial point on the characteristic $(z_1,t_1)$. Thus, I want to solve for the characteristic starting at $(z_1,t_1)$, such that the evolution of $F$ along the characteristic is described by $dF/ds$. This characteristic is illustrated in Fig.~\ref{fig:mol} where it traces a path in the space defined by $z$ and $t$. The chain rule allows me to rewrite the left hand side of Eq.~(\ref{eqn:char_wave}) as:
\begin{equation}
  	\frac{dF}{ds}=\frac{dz}{ds}\frac{\partial F}{\partial z}+\frac{dt}{ds}\frac{\partial F}{\partial t}.
\end{equation}
The equations have already been scaled such that the characteristics have a simple form. If I set
\begin{equation}\label{eqn:char_coefs}
	\frac{dz}{ds}=1 \quad \mathrm{and} \quad \frac{dt}{ds}=1,
\end{equation}
then \ref{eqn:char_wave} gives:
\begin{equation}\label{eqn:char_ode}
  \frac{dF}{ds}=(\mathscr{L}+\mathscr{N})F.
\end{equation}

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[scale=1]{Figures/mol.pdf}
  \end{center}
  \caption[Method of characteristics]{The characteristic for Eq.~(\ref{eqn:char_wave}), illustrated in $z$-$t$ space, can be described by $dF/ds$. Using variables with physical units, $\Delta z=(c/n)\Delta t$, hence the characteristic is the light-line for a propagating wave.}
  \label{fig:mol}
\end{figure}

Equation~(\ref{eqn:char_coefs}) gives $z=s+z_1$ and $t=s+t_1$. Solving these for $s$ and equating yields:
\begin{equation}
  z-z_1=t-t_1,
\end{equation}
or simply $z=t$ where I have taken $(z_1,t_1)=(0,0)$.

These characteristics have a simple physical interpretation, by momentarily reverting to variables with physical units via Table~\ref{table:normalized}, we see that 
\begin{align}
  \frac{z}{L}&=\frac{t}{t_r}=\frac{tc}{n_0 L}\\
  z&=\frac{c}{n_0}t,
\end{align}
also illustrated in Fig.~\ref{fig:mol}. Hence, the characteristics are the lines through space-time that describe a wave propagating along the $z$ axis with phase velocity $\pm c/n_0$. In the case of wave equations, the characteristics always have this interpretation, and thus the method of characteristics can generally be used to reduce PDEs describing wave propagation to ODEs along the light-lines of each propagating wave.

Numerical integration may thus take place along the characteristics of the wave equation by taking each numerical step in time and space simultaneously. In the case presented here, this allows the simplification of Eq.~(\ref{eqn:model_eqns_F}) as
\begin{equation}
  \label{eqn:df_dt}
  \frac{dF}{dt}=(\mathscr{L}+\mathscr{N})F,
\end{equation}
where $z$ is replaced with $t$ corresponding to integration along the characteristic $z=t$. The following section describes how numerical integration of Eq.~(\ref{eqn:df_dt}) is achieved via the split-step beam propagation method.

% section Numerical Methods (end)

\section{Split-step beam propagation method} % (fold)
\label{sec:split_step_beam_propagation}

The method of characteristics has reduced the partial differential equation given by Eq.~(\ref{eqn:char_wave}) to an ordinary differential equation. The primary difficulty in solving numerically this set of equations is that they contain both linear and nonlinear terms. The split-step BPM \cite{Fleck_1976aa} is a numerical technique for efficiently integrating the linear and nonlinear portions of a wave equation by alternately treating linear diffraction in the Fourier domain and the nonlinear terms in the spatial domain.

The foundation of the BPM is that the Maxwell equations for counter-propagating beams in a nonlinear medium can be approximated by a set of $m$ thin slices separated by regions of empty space.\footnote{I have chosen to denote the number of slices $m$ in order to avoid any confusion with the index of refraction. In the literature, $n$ is commonly used for both.} Waves propagating through this sequence of spaces and slices experience only linear diffraction in the empty space and only a nonlinear phase shift within each slice. This model medium is illustrated in Fig.~\ref{fig:app_slices}.

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[scale=1]{Figures/bpm.pdf}
  \end{center}
  \caption[A model medium for the split-step BPM]{A model medium of length $L$ illustrating the split-step BPM. Slices of nonlinear media (solid black lines) separated by free-space (grey) of length $\Delta z$ combine to approximate the effects of a medium that exhibits both linear diffraction and nonlinear wave mixing. The dashed grey lines indicate the boundary between each unit cell consisting of of free space, nonlinear slice, and free space.}
  \label{fig:app_slices}
\end{figure}

The approximation of the medium as a series of slices originates from approximating the solution to Eq.~(\ref{eqn:df_dt}) as \cite{Penman_1990aa,Agrawal_2001aa}
\begin{equation}
  \label{eqn:bpm}
  F(x,z_1,t_1)=\exp\left(\mathscr{L}\frac{\Delta t}{2}\right)\exp\left(\mathscr{N} \Delta t\right)\exp\left(\mathscr{L}\frac{\Delta t}{2}\right)F(x,z_0,t_0)+\mathcal{O}(\Delta t^3),
\end{equation}
where $\Delta t=t_1 - t_0$, and $z_0=z(t_0)$, with a similar equation for $B$.

It is straightforward to apply the nonlinear operator in the spatial domain because the operator contains only multiplication operations. From Eq.~(\ref{eqn:nonlinear_op}), a single time-space step through a nonlinear slice is calculated via
\begin{equation}
  F(x,y,z_1,t_1)=\exp(\mathscr{N}L_m)F(x,y,z_0,t_0)=\exp[i(|F|^2+2|B|^2)L_m]F(x,y,z_0,t_0),
\end{equation}
where $L_m=L/m$ for $m$ slices.

The linear operator could be applied in the spatial domain; however, such an operation requires calculating a second-order spatial derivative, which is computationally intensive. Fortunately, the linear operator can be applied in the Fourier domain very efficiently via the Fast-Fourier Transform algorithm \cite{Fleck_1976aa}. For the Fourier transform $\widetilde{F}(\kappa_x,\kappa_y,z,t)$ of $F(x,y,z,t)$, Eq.~(\ref{eqn:linear_op}) becomes
\begin{equation}
  \frac{d\widetilde{F}}{dz}=\frac{i}{4\pi\mathscr{F}}(\kappa_x^2+\kappa_y^2)\widetilde{F},
\end{equation}
with solution
\begin{align}
 \widetilde{F}(K_x,K_y,z_1,t_1)&=\exp(-i\theta_n)\widetilde{F}(\kappa_x,\kappa_y,z_0,t_0)\label{eqn:linear_operator}\\ \theta_n&=\frac{(\kappa_x^2+\kappa_y^2)\Delta z}{4\pi\mathscr{F}}\label{eqn:theta_n},
\end{align}
where $\Delta z=L/(m-1)$, with $m$ slices.

A subtlety of the solution described in Eq.~(\ref{eqn:bpm}) is that the approximation is only exact for a physical system where the length of the nonlinear medium is not equal to the total length of propagation. This is illustrated in Fig.~\ref{fig:app_slices} where there is an extra half-step of free-space propagation on either side of the nonlinear slices. The nonlinear medium I use in my experiments is contained within a glass cell of length $L=5$~cm, so including the extra half-step on each side of the simulated medium is equivalent to measuring the field a small distance away from the cell. However, the measurements of the switch response I make experimentally are all in the far-field; hence, this approximation is suitable for simulating my experimental system. If I were to simulate the effects of placing this medium within a cavity, extra consideration must be paid to the effect of additional regions of free-space propagation.

Another consequence of approximating a continuous medium with a set of discrete slices is an artificial instability that can arise due to the inherent periodicity of the slices. In the thin-slice model, each beam experiences phase modulation within a nonlinear slice and the free-space propagation can then convert this phase modulation into an amplitude modulation via diffraction. Penman \etal \cite{Penman_1990aa} show, by linear stability analysis, that specific values of the transverse wavenumber $K$ experience no amplitude modulation (\emph{i.e.}, they are only phase-shifted) due to the specific periodicity of the model system. This can result in an instability for large-$K$ that does not correspond to an instability in the continuous model. To avoid this artificial instability, spatial filtering is implemented that excludes transverse wavevectors that satisfy $K>(\pi k/\Delta z)^{1/2}$, or using the notation of Eq.~(\ref{eqn:theta_n}), $\theta_n>\pi/2$. During a discussion of this issue, Dr. William Firth clarified for me that this filtering was implemented by Chang \etal \cite{Chang_1992aa}. This filtering has been applied in my simulations, and I do not observe the large-$K$ instability.

A final step in modeling hexagon formation numerically is a source of symmetry breaking required to overcome the square symmetry of the numerical grid. The model implemented by Chang \etal \cite{Chang_1992aa} successfully demonstrated two sources of symmetry-breaking, noise added to the initial gaussian beams, and anisotropy in the numerical grid. I have observed hexagon formation using either of these techniques in addition to a third technique where I inject a weak off-axis probe beam for a few transit times very early in the course of the simulation. This weak beam serves to seed the instability and also selects the initial pattern orientation.

% section split_step_beam_propagation (end)

\section{Numerical Results} % (fold)
\label{sec:numerical_results}

The previous work of Chang \etal shows that three-dimensional simulations of counterpropagating gaussian beams exhibit transverse patterns. Specifically, they demonstrate that generalizing the model of Firth and Par\'e to two transverse dimensions leads to the formation of hexagonal patterns in agreement with a large number of experimental results (see \cite{Chang_1992aa} and References therin). Chang \etal find the formation of hexagons for pump powers just above the threshold predicted by Firth and Par\'e and well below the threshold predicted by Yariv and Pepper. 

It must be made clear that neither the simulations conducted by Chang \etal nor the simulations I conduct, are intended for direct comparison with experimental results. This would require refinement of the model for the nonlinear response of the medium to accurately describe the nonlinear response of rubidium vapor (see Chapter~\ref{cha:patterns} for a discussion of the limitations of the Kerr model).

\subsection{Pattern Formation} % (fold)
\label{sub:pattern_formation}

The primary result presented by Chang \etal is the formation of hexagonal patterns in a three-dimensional model of gaussian beams counterpropagating through a medium exhibiting Kerr nonlinearity. Their simulations are conducted with $\mathscr{F}=63.7$ and $IL=0.565$ where the threshold for plane-wave pattern formation predicted by Firth and Par\'e is $IL\simeq0.45$. Therefore, Chang \etal simulate pattern formation for pump beams that are 25\% above the minimum plane-wave threshold.

I have conducted simulations with a wide range of $\mathscr{F}$ between 63.7 and 4, where my experimental conditions correspond to $\mathscr{F}=5.3$. Simulations in this range all exhibit hexagonal pattern formation and reproduce the results of Chang \emph{et al.}; I use this as one criterion for confirming my computational implementation of the BPM.\footnote{Other criteria for confirming that I have correctly implemented the BPM are simulating gaussian beam propagation \cite{Siegman_1986aa}, and nonlinear self-focusing \cite{Banerjee_2000aa}, both of which are correctly reproduced using the appropriate parameters.} The primary difference between the results for different $\mathscr{F}$ is simply the opening angle $\theta$ of the cone describing the generated light. From the analysis of Firth and Par\'e, the instability threshold is lowest for $K^2L/2k\simeq3$. Given $\mathscr{F}=w_0^2/\lambda L$, and $\kappa=w_0K$, the minimum instability threshold corresponds to a dimensionless transverse wavevector of $\kappa^2=12\pi\mathscr{F}$. Clearly larger $\mathscr{F}$ results in a larger transverse wavevector generated by the instability; hence, my observation that small $\mathscr{F}$ results in small $\theta$. In order to simulate the geometry of my experiment, the results presented below are of simulations where $\mathscr{F}=5.3$ and $IL=0.565$ ($\sim25$\% above threshold).

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[scale=1]{Figures/sim_frames.pdf}
  \end{center}
  \caption[Patterns and switching shown by numerical model]{Numerical simulation of counterpropagating gaussian beams shows ring and hexagon pattern formation in the far field.}
  \label{fig:sim_frames}
\end{figure}

Images of the far field pattern generated in a typical run of my simulation are shown in Fig.~\ref{fig:sim_frames} where the time corresponding to each frame is indicated in units of the transit time $t_r=n L/c$. In the initial frame of Fig.~\ref{fig:sim_frames}, the transmitted forward pump-beam is visible at the center, and the weak off-axis perturbation is visible to the right. This perturbation is used to quickly induce hexagonal pattern formation. Without the initial perturbation, hexagons are observed after 100-150 transit times. At $t_r=17$, the field that is conjugate to the perturbation, and due to forward four-wave mixing, is visible to the left of the central pump. The Fast-Fourier Transform (FFT) is used to compute these far-field images, and one artifact of this is a large-narrow peak corresponding to zero spatial frequency (DC). This feature is located in the center of the frame and is much narrower than the transmitted pump beam. The dark dot in the central spot of the first two frames is the result of numerical filtering that I use to remove this component of the image. At $t_r=23$, a ring pattern has formed that is replaced by hexagons at $t_r=53$. The perturbing beam is turned off at $t_r=35$ and is not visible at $t_r=53$. It is interesting to note that the ring pattern, predicted by naively generalizing the models of Yariv and Pepper or Firth and Par\'e to cylindrically symmetrical transverse dimensions is a transient solution that appears early in the development of the off-axis patterns. The ring is not a stable solution for the system in the presence of symmetry breaking, due in this case to the initial perturbation beam, and the ring breaks up into six spots after a short time. The second row of frames shown in Fig.~\ref{fig:sim_frames} were collected after the application of an off-axis switch-beam, which turns on at $t_r=85$, and are discussed in the next section.

% For the geometry of my experiments, $L=5$~cm, this timescale is $t_r\simeq0.16$~ns, yet as I have shown, the typical timescale exhibited by my experiment is much larger, on the order of $\mu$s. This is one indication that the Kerr model will not quantitatively represent my experimental system, however, the patterns I observe, and more importantly, the effect of an external perturbation on those patterns is qualitatively similar.

% subsection pattern_formation (end)

\subsection{Switch Response} % (fold)
\label{sub:switch_response}

In my simulations, much like in my experiments, I observe that injecting a weak switch beam into the nonlinear medium after hexagons have formed causes the hexagon pattern to rotate such that a bright spot is aligned to the direction of the switch beam. This rotation is illustrated in the lower four frames of Fig.~\ref{fig:sim_frames}. The switch beam is applied at $t_r=85$, and becomes visible between the two right-side spots at $t_r=101$. For the frames shown, the switch beam power is $P_s=10^{-4}P_p$, where $P_p$ is the power of each of the counterpropagating pump beams. At $t_r=150$, the counterclockwise rotation of the pattern can be observed and continues until the end of the simulation at $t_r=300$ where the pattern has rotated such that the locations that were previously bright are now dark.

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[scale=1]{Figures/sim_apertures.pdf}
  \end{center}
  \caption[Location of apertures used in numerical models]{The location of the on- and off-state apertures are indicated relative to the initial hexagon pattern that forms at $t_r=53$. The on-state aperture (upper square) is located opposite the applied switch beam, and the off-state aperture (lower square) transmits the spot immediately counter-clockwise from the on-state aperture.}
  \label{fig:sim_apertures}
\end{figure}

As in the switch I demonstrate experimentally, the patterns generated in this simulation can be spatially filtered in order to define two or more output channels. Figure~\ref{fig:sim_apertures} indicates the location of the apertures used to filter the numerical results. Square apertures are used for numerical efficiency, but the results are not expected to differ if they are replaced with round apertures. The power transmitted by these apertures is calculated by summing the simulated intensity values within each aperture. For four simulation runs, each with different switch-beam power, the power transmitted through the on- and off-state aperture as a function of transit time $t_r$ is shown in Fig.~\ref{fig:sim_response}(a) and (b), respectively.

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[scale=1]{Figures/sim_response.pdf}
  \end{center}
  \caption[The power transmitted by simulated apertures.]{The power transmitted by apertures in the numerical model exhibits switching behavior that is similar to the experimental system. The response of the on- and off-state aperture transmission is shown for four levels of switch-beam power. The switch beam is turned on at $t_r=85$, indicated by the arrow in (a). As the switch-beam power decreases, the simulation exhibits slower response, \emph{i.e.} slower pattern rotation. The switch-beam power (in units of pump-beam power) corresponding to these four traces are $10^{-4}$ (solid black), $2.5\times10^{-5}$ (lg. dash blue), $4\times10^{-6}$ (sm. dash red), and $1\times10^{-6}$ (dash-dot green). The horizontal dotted line indicates the threshold used to calculate response times for the simulated switch.}
  \label{fig:sim_response}
\end{figure}

After the initial transients in the pattern formation, the power in the off- and on-spots stabilize within 50 transit times. At $t_r=85$, the switch-beam is applied and the pattern begins to rotate, transferring power from the off-state aperture to the on-state aperture. For $P_s=10^{-4}P_p$, complete rotation occurs within 200 transit times. For lower switch beam power, the pattern rotates more slowly as the remaining traces show in Fig.~\ref{fig:sim_response}. To compare the change in response time observed in the simulation to that observed experimentally, I measure the response time of the simulated switch as the time between the application of the switch beam ($t_r=85$) and the threshold crossing for the on-spot. The threshold, also shown in Fig.~\ref{fig:sim_response}, is chosen to roughly correspond to the level of the experimental threshold.

\begin{figure}[htbp]
  \begin{center}
    \includegraphics[scale=1]{Figures/sim_response_time.pdf}
  \end{center}
  \caption[Simulation of the switch response time.]{Simulation of the switch exhibits an increase in response time for decreasing power that is qualitatively similar to experimental observations. To facilitate comparison to Fig.~\ref{fig:photon_number}(a), the horizontal axis has high switch-beam power to the left and low switch-beam power to the right.}
  \label{fig:sim_response_time}
\end{figure}

The response time of the simulated switch data shown in Fig.~\ref{fig:sim_response_time} ranges from 40 transit times to 210 transit times, as shown in Fig.~\ref{fig:sim_response_time}. For comparison, the transit time of the 5~cm vapor cell used in my experiment is 160~ps, so the simulated response times would correspond to experimental values of 6.4~ns and 33.6~ns respectively. As described in Chapter~\ref{cha:switch}, I observe response times between 2 and 4~$\mu$s experimentally, so it is clear that the numerical model does not agree quantitatively with my experimental observations. However, the simulated response time does exhibit a sharp increase in the limit of low switch-beam power, which is qualitatively similar to my experimental observations. This increase in response time for weak inputs indicates that the switch undergoes critical slowing down, which suggests that the orientation of the pattern exhibits multi-stability between the preferred orientations.

% In order to increase the quantitative accuracy of the model it would be necessary to include the nonlinearity due to optical pumping effects. Optical pumping in Rb typically occurs with a timescale on the order of $\mu$s, and thus correspond well to my experimental observations of the switch response time.

Another notable feature of these numerical results is that, despite the limitations of the model, the amount of switch-beam power, relative to the total pump power, required to rotate the pattern is of the same order of magnitude as what I observe experimentally. For reference, Table.~\ref{tab:exp_sim} shows the correspondence between the normalized switch-beam power $P_s/P_p$ used in the simulations presented above and the experimental values, based on total pump-beam power, from my experiments, of $P_p=560~\mu$W. As an example, the third curve in Fig.~\ref{fig:sim_response} corresponds numerically to $P_s=4\times10^{-6} P_p$. In my experiment, this switch-beam to pump-beam power ratio would imply a switch-beam power of 1~nW. In my experiment, the switch typically operates between 1~nW and 50~pW. Therefore, the sensitivity exhibited by my experiment is largely described by this model.

\begin{table}
  \begin{center}
  \begin{tabular}{c|c}
    $P_s/P_p$&$P_s$[nW]\\
  	\hline
  	$1\times10^{-4}$&26\\
  	$2.5\times10^{-5}$&6.5\\
  	$4\times10^{-6}$&1\\
  	$1\times10^{-6}$&.26\\
  \end{tabular}
  \end{center}
  \caption[Correspondence between $P_s/P_p$ and $P_s$ in nW]{The correspondence between $P_s/P_p$ and $P_s$ in nW based on 560~$\mu$W of total pump power.}
  \label{tab:exp_sim}
\end{table}

There are certainly features of the experiment that these simulations do not capture. In the first case, absorption is neglected in assuming a Kerr-type nonlinearity. One consequence of this is that misalignment of the pump beams in the simulation does not serve to reduce the number of spots generated. This is likely due to the fact that, without absorption, there is no loss experienced by the less-favored hexagonal components and even for large pump-beam misalignment, the pattern remains a hexagon. Simulations that include misalignment of the forward pump beam exhibit hexagonal pattern formation in addition to fluctuations in the pattern orientation and a near-field pattern flow \cite{Honda_1995aa}. Just as for well-aligned pump beams, the switch beam also causes pattern rotation when the forward pump beam is misaligned, and the switch response time diverges near zero switch-beam power in the misaligned case as well.

As I discuss in the following Chapter, symmetry breaking may be responsible for pinning the orientation of the pattern. However, including pump-beam misalignment in these simulations does not have this effect: even weak switch beams cause pattern rotation and there does not seem to be a critical switch-beam power required to rotate the pattern. Refinement of the model to include absorption and saturation may improve the agreement between my experiment and the simulation. Furthermore, because I have assumed a medium with an instantaneous nonlinear response, the only timescale in the Kerr model is the transit time. This leads to significantly faster switch response in the model compared to my experimental observations. To quantitatively model the response time requires a more refined model of the nonlinear interaction that includes optical pumping effects and the associated time scales.

% subsection switch_response (end)

% section numerical_results (end)

\section{Summary} % (fold)
\label{sec:sim_summary}

Pattern formation in nonlinear optical systems has been well-studied via experiment, theory and numerical modeling. Previous simulations, based on the Kerr model described by Firth and Par\'e, have been conducted by Chang \etal and exhibit the formation of hexagon patterns for a wide range of simulation parameters. I extend this result by showing that not only are hexagonal patterns generated by this model, but the model also describes pattern rotation induced by an external switch beam. Furthermore, this model successfully reproduces the general relationship between the switch response time and the switch-beam power that I observe experimentally.

The only timescale included in the model is the transit time, so the fact that the response time depends on the strength of the perturbation may be a universal feature of all-optical switching with transverse patterns.

The qualitative agreement between this simple model and my experimental observations indicate that the transparent Kerr medium, although different from rubidium vapor in important ways, describes many of the features that make transverse optical patterns useful to applications in all-optical switching.

% section sim_summary (end)

